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This is a summary report of a two-year grant. We started to develop two- and three-dimen- 
sional MHD equilibrium models for Earth's magnetosphere. Our original proposal was moti- 
vated by realizing that global, purely data-based models of Earth’s magnetosphere are inade- 
quate for investigating the underlying plasma-physical principles according to which the magne- 
tosphere evolves on the quasi-static convection time scale (see Section I). 


During the two years of grant NAGW-1781, we established complex numerical grid-gener- 
ation schemes for a three-dimensional Poisson solver, and we succeeded in coding a rather ro- 
bust Grad-Shafranov solver for high-beta MHD equilibria. Thus, we were able to calculate the 
effects of both the magnetopause geometry and boundary conditions on the magnetotail current 
distribution. That work is summarized in Section II. 
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I. Objectives: 

Magnetospheric MHD Equilibrium Configurations: 

The problem of plasma equilibrium for Earth's magnetosphere is of fundamental importance. 
Any attempt to model processes on the quasi-static time scale, such as magnetospheric con- 
vection, or the interaction between the ionosphere and the magnetosphere in terms of field 
aligned Birkeland currents, requires that the magnetic field model satisfy the condition of quasi- 
static MHD equilibrium. 


Conventional magnetospheric B-field models which only satisfy Maxwell’s equations merely 
reproduce the observed magnetic field structure in terms of empirically determined current dis- 
tributions (magnetopause currents, tail currents, ring current). Those models do not include the 
underlying plasma physics that would explain any particular magnetic field configuration. 

In the slow-flow (quasi-static) MHD approximation, the convection flow field is decoupled 
from the the force balance equation; thus one has to solve the equilibrium equations 


j x B = V • P 

(1) 

V x B = |io j 

(2) 

V • B = 0 

(3) 


with appropriate magnetopause boundary conditions for the magnetic field, and plasma sheet 
boundary conditions for the thermal plasma pressure. 

Ideally, one would like to have a three-dimensional magnetospheric equilibrium model that 
would allow us to "play around" with various magnetopause boundary conditions, boundary 
shapes, and thermodynamic conditions for both isotropic and anisotropic thermal plasma. Such 
an ideal model would then allow us to see whether the system reaches a steady state or con- 
serves one or the other thermodynamic quantity in the course of its time evolution. Also, MHD 
equilibrium configurations can be tested for stability, a problem which is intimately related to the 
problem of magnetospheric substorms! 

Thus we began two years ago with a rather ambitious (NASA supported) research program 
that eventually will provide us with an "ideal" MHD equilibrium model explained above. Nu- 
merical equilibrium codes adopted from fusion plasma physics cannot simply be used in a high 
plasma beta ((3 » 1) environment such as Earth's magnetotail. We therefore started from the 
very beginning from building codes that are based on sophisticated numerical grid generation 
schemes (see Section II. 1). Following that, we were able to develop a three-dimensional Pois- 
son solver for Neumann boundary value problems (Section n.2), and a Grad-Shafranov solver 
for two-dimensional equilibria (Section II.3). 
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Those new codes have very recently allowed us to address interesting magnetospheric 
physics questions regarding the influence of the magnetopause shape and boundary conditions 
on the self-consistent magnetotail currents (Section II.4). 


II. Accomplishments from Previous Grant NAGW-1781: 

We spent most of the our budgeted time during the last two years with the development of 
rather sophisticated numerical codes we describe in this section. Since the material discussed 
below has not been published yet, we explain our main accomplishments in same detail. 

1. Numerical Grid Generation: 

During the first funding cycle, Frank Toffoletto has implemented a computational grid- gener- 
ation scheme [Thompson et al., 1985] that allows us to formulate numerical difference equations 
such that derivatives can be calculated with high accuracy. Numerical grid methods have the ad- 
vantage that the computational domain becomes a uniformly spaced rectangular grid that greatly 
simplifies the implementation of the the boundary conditions: the grid boundary exactly corre- 
sponds to the physical boundary of the system. 

The proper choice of an optimal non-orthogonal, non-cartesian grid is important for the con- 
vergence behavior of an MHD equilibrium scheme. The optimal choice of the grid often depends 
on the problem at hand. However, maximum orthogonality in the numerical grid pattern is desir- 
able in order to minimize truncation errors in the difference equations. 

We decided to use a grid generation algorithm of maximum flexibility: Our magnetospheric 
grids were computed by the algebraic grid generation program TBGG. TBGG was developed at 
the Los Alamos National Laboratory; it allows interactive user control of the grid structure via 
modification of such features as grid spacing and orthogonality. The program generates compu- 
tational grids of rectangular structure. Such grids are convenient in "finite difference" or "finite 
volume" solutions of boundary and initial value problems. Computational grids of rectangular 
structure are topologically equivalent to a square: for each grid point in the physical space (x, y, 
z coordinates) there is a corresponding point in rectangular computational space (Q, £), ^). 


2. The Three-Dimensional Poisson Solver: 

Based on the grid generation method described above, we were able to establish a code that 
solves Poisson's equation (or Laplace's equation as a special case) in three dimensions. The 
conservative version of the Laplacian operator is given by 

1 3 3 

v 2 ( P=~rSS[a 1 -(a J( P)^H i ( 4 ) 

g i=i j=i 

where ()j^ is the differentiation with respect to the coordinate t, 1 . By using the metric identity 
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X ( g 1/2 a ')^=0 


i=l 


( 5 ) 


the Laplacian (4) reduces to the form 


3 3 

V 2 cp=-4-SI(g 1/2 g ij cp^ (6) 

g ' i=i j=i 

which is conservative in the Jp direction. The term g 1/2 is the Jacobian defined from the volume 
element dV = g 1/2 d^p dip d^ k , and g ij is defined as the dot product 

g ij = V^-V^ (7) 


The expression Vlp is determined from 


1 5r dr 


(i, j, k cyclic) 


where r is the position vector of the grid in physical space at the location (Ip, Ip, 
The magnetic field components are determined from the conservative expression 


B = - V(p 


, 1/2 t 


S (g 1/2 a i 9 )$' 


i=l 


( 8 ) 


( 9 ) 


Grid Discretization : 

To solve the Poisson equation on a numerical grid, we tried two methods of discretization, 
namely, the "finite difference", and the "finite volume" approach. 

In the "finite difference" approach, the scalar potential <j> and its derivatives are calculated 
at the grid locations. In two dimensions, we have used both an average of a forward-backward 
scheme and a second order scheme [Thompson et al., 1985] for solving Laplace’s equation. A 
three-dimensional generalization of those two schemes is straightforward. 

Both schemes are symmetric with respect to the grid points which ensures hat the second 
derivatives satisfy the relation 


a 2 cp _ a 2 tp 


( 10 ) 


In the "finite volume" approach, the scalar potential $ and its derivatives are calculated at 
the center of a cell which is defined as a square surrounding the grid point (a cube in three di- 
mensions); the derivatives of (J) are calculated at the boundaries between the cells. The "finite 



5 


volume" scheme also satisfies (10) and reduces to the well known second order scheme for a 
purely rectangular grid (X = 1), namely, 

V 2 <p = (p i+ X,j +(pi-Xj +cpi, j+ x+(pij.x-4 <p»j (11) 


Magnetic Boundary Conditions: 

Magnetospheric physics problems typically require the solution of Neumann boundary value 
problems at the magnetopause [e.g., Voigt, 1981]. For Neumann boundary conditions, the nor- 
mal component of the magnetic field is given by 


B„ = .-^L*B (12a) 

I V4 J I 

By using (5) and (9), equation (12a) reduces to 

, 3 

Bn = -^=Eg' JA 4i d2b) 

Vg-® i=l 

A "finite volume" approach is more convenient for the solution of Neumann boundary value 
problems. The implementation, however, is a little more complex, because it requires the calcu- 
lation of the metric coefficients at the cell faces as well as the cell centers. For a (m x n) grid 
there are (n-1) x (m-1) cells and (n+m-1) ‘ghost’ cells that are needed to impose the Neu- 
mann boundary conditions. These ‘ghost’ cells lie outside the range of the computational box, 
so that the normal component at the magnetopause, for example, is given by 

1 1 12 

Bn = rj, ( ijmax-l/2((Pi+l jmax + 

g ij-1/2 4 

(Pi+1 jmax-1 ~ (Pi-1 jmax _ (pi-1 jmax-l) + G i jmax-l/2((Pijmax — *Pi jmax-l)) (13) 

which provides a second-order accurate, discretization of the magnetic normal component at the 
magnetopause. 

Solution of Poisson’s Equation: 

We briefly mention that the algorithms sketched above lead to the solution of the matrix 
equation 


MA = S (14a) 

A = M _1 S 


(14b) 
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For a two-dimensional problem on a (n x m) grid, M is a sparse (n x m) 2 matrix, A is a 
vector of order (n x m) corresponding to the potential in equation (4), and S is a vector denot- 
ing the source term for (4) as well as the magnetic boundary conditions. 

In three dimensions, owing to the large sizes of the matrices, an iterative scheme must be 
employed to evaluate (14b). We have utilized the software package PCGPAK available on the 
CRAY-YMP at NCSA for that purpose. The package uses a conjugate gradient method 
adopted for non-symmetric matrices. The successful inversion of the matrix M depends on the 
choice of grid which effects the size and number of the off-diagonal terms (i.e., terms involving 
gy , where i*j). A second-order difference scheme in three dimensions typically involves over 
30 off-diagonal terms. The ultimate solution is usually obtained after several trial grids. 

The two-dimensional version of our codes produce sufficiently small matrices that direct ma- 
trix inversion techniques can be employed, for this we use the matrix package SMPAK which is 
based on the Yale Sparse Package and is also available on the CRAY-YMP at NCSA. Again, 
care must be taken in choosing the appropriate grid. 


3. The Two-dimensional Grad-Shafranov Solver: 

Two-dimensional magnetospheric MHD equilibria for isotropic plasma pressure can be cal- 
culated by solving the Grad-Shafranov equation, 

o dP 

V 2 a + po = — 1*0 j dipole (15) 

da 

This equation is equivalent with the equilibrium equations (1) to (3) for V • P — > VP. The 
function a(x, z) is a (scalar) magnetic flux function and satisfies the conditions B • Voc = 0 
and V • B = 0. The variables x and z are used in a (GSM) coordinate system where the x 
axis points toward the sun, and the x - z plane contains the planetary dipole. In the isotropic 
limit, the plasma pressure P(a) is constant on magnetic field lines. The plasma currents are 
given by 


j y «X) = — (16) 

da 

and the magnetization current which describes the planetary magnetic dipole depends on the 
dipole tilt angle \|/ according to 


P0 jdipole = - Md [^-S(x) 5(z) cos \|/ - 5(x) ^-5(z) sin \|/ ] (17) 

where Md is the magnetic dipole moment. The dipole is located at the origin of the (GSM) 
coordinate system. 
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Analytic solutions of (15) which include the dipole term (17) exist only for a very specific 
pressure function, namely, 


P(CC) = — !— k 2 a 2 (18) 

2 M-0 

and boundary conditions specified on a rectangular magnetopause [e.g., Voigt, 1986]. Thus we 
decided to develop a so-called Grad-Shafranov solver for high beta-plasma, i.e., a numerical 
scheme that solves the equilibrium equation (15) for general pressure functions P(a), and for 
realistic magnetopause geometries and boundary conditions (more details follow in Subsections 
II.4.1 and H.4.2). 

For solving the equilibrium equation (15) it is convenient to split the magnetic flux function 
into two separate components, 


CX — Ot^ipoie + Otj (19) 

where a dipo i e is the flux function for the Earth's dipole field and (Xj is the flux function for the 
plasma currents. The flux function for the dipole satisfies the relation 


V 2 


a dipole — M-Oj dipole 


Thus equation (15) reduces to 


V 2 cxj + = 0 

J da 


( 20 ) 


( 21 ) 


If we use the pressure function (18), for example, then equation (21) reads 

V 2 aj + jj-o k 2 exj = - it ok 2 a dipole (22) 

The solution of (22) involves solving a matrix inversion problem of type (14), where S is the 
right side of (22). At the far-tail boundary, we have assumed for simplicity that 

= 0 (23) 

K 1 

The "closed" magnetosphere, i.e., B n = 0 at the magnetopause, requires Dirichlet bound- 
ary conditions which are specified by the simple condition aj = - a dipole . 

More interesting is a generalized "open" magnetosphere. General non-zero boundary con- 
ditions can be formulated in terms of 


n • B — B n — /(x m p) 


(24a) 
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where (x mp ) denotes locations at the magnetopause. Using the relations (5), (9), and (12a), 
and noting that B = Vaxy we obtain the magnetic normal component (24a) at the magne- 
topause in the form 


B n = - 


9otj 

W 

. 22\U2 

(g g ) 


The integration of (24b) gives the boundary values for the magnetic flux function, 


= - OC djpole - f B n (g g 22 ) 


22 ' m dV 


(24b) 


(25) 


where denotes the location where a = 0. is chosen to be the subsolar point. Following 
the theory of Toffoletto and Hill [1989] for the "open" magnetosphere, the magnetic normal com- 
ponent, B n = /(x mp ), can be specified by the relation 

Vx(v t x B n ) = 0 (26) 

which is subject to an input velocity function 

v t = v sw — + v a (l-— ) (27) 

Z oo Zoo 


where v sw is the solar-wind velocity, v a is the Alfvdn velocity, and is the far-tail radius. 
Equation (27) approximates the tangential magnetosheath flow away from the subsolar point 
plus an additional merging outflow which is of the order of the Alfven speed. In two-dimensions, 
(26) reduces to the expression 


v t B n = v sw B n oo = constant 

B noo is the normal component in the far-tail. Then B n becomes 

B noo 


(28) 


B n (z) = 


Z oo V sw Z oo 


(29) 


The magnitude of the constant B noo can be estimated from the polar cap co-latitude (p 0 , the tail 
length L, and dipole moment M as 


B 


n °° 


2 M sin 90 


L 


(30) 
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because a constant (<po B n ) is much smaller than the corresponding normal component magni- 
tude for a three-dimensional magnetosphere. The corresponding expression for (30) in three 
dimensions is given by 


7t M sin 2 cpo 
2 L R 


where R is the tail radius. The ratio of the normal component is then 


2-D , t> 3-D 7C Sin (po 


n©o / ~ 


4 R 


(31) 


(32) 


which is much smaller than unity. 

4. A Bouquet of Physical Results: 

With the Grad-Shafranov solver in place, we were able to address a number of interesting 
physical questions regarding the influence of the plasma pressure and of magnetopause bound- 
ary conditions on the equilibrium magnetotail structure. 

4.1 The Closed Magnetosphere: 

The "closed" magnetosphere, /(x mp ) = 0, remains as a special case in a class of more gen- 
eral solutions to the boundary conditions (24a). That simple boundary condition allows us to 
demonstrate quite nicely the influence of the plasma pressure and the magnetopause shape on 
the magnetotail configuration. 

I n fl uence of the magnetoml plasma pressure: 

In a feasibility study for this proposal, we chose for the plasma pressure the functional form 
(18). The parameter k in (18) is a measure of the plasma pressure. As the thermal pressure 
increases (k increases), owing to an unspecified plasma supply mechanism, tail currents de- 
velop and stretch the magnetotail field lines accordingly. Of course, the pressure function (18) 
can hardly be justified from a thermodynamic point of view. Thus we began trying more realistic 
pressure functions which will allow us to investigate series of magnetospheric equilibrium 
states that are compatible with magnetotail convection. This will then allow us to study the un- 
derlying thermodynamics of the magnetosphere. 

Influences/ t he m agnetopause sh ap. ei 

We found that the magnetopause geometry strongly influences the overall magnetotail struc- 
ture; the field lines change from a more dipolar configuration to a very stretched tail configura- 
tion without changing the plasma pressure function P(a). We were able to demonstrate that 
the solar wind pressure itself, via its influence on the magnetopause geometry, affects the very 
structure of the magnetotail plasma currents. This important effect can be investigated on the 
level of the MHD equilibrium theory; it does not reveal itself in traditional empirical models that 
only satisfy Maxwell's equations. 
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4.2 The Open Magnetosphere: 

It is well known that many aspects of magnetospheric dynamics have to be modeled by 
means of an "open" magnetosphere, i.e., for general magnetic boundary conditions /(x m p) in 
(24a). From the very beginning, we therefore coded the new Grad-Shafranov solver (see Sec- 
tion II.3) such that general boundary conditions can be included. 

Influence of the magnetopause boundary conditions: 

The "open" magnetosphere possesses a small southward B z component along the tail 
magnetopause. Such a B z component was calculated according to the theory of Toffoletto and 
Hill [1989] we sketched briefly in equations (24) to (32). We were able to demonstrate that a 
southward turning of the IMF can, indeed, stretch the magnetotail and increase the tail current 
density accordingly. This effect is extremely important for the problem of substorm dynamics. 


References: 

Thompson, J. F., Z. U. A. Warsi, and C. W. Mastin, Numerical Grid Generation, Foundations 
and Applications, North-Holland, Amsterdam, 1985. 

Toffoletto, F. R., and T. W. Hill, Mapping of the solar wind electric field to the Earth's polar caps, 
J. Geophys. Res., 94, 329-347, 1989. 

Voigt, G.-H., A mathematical magnetospheric field model with independent physical parameters, 
Planet. Space Sci., 29, 1-20, 1981. 

Voigt, G.-H., The shape and position of the plasma sheet in Earth's magnetosphere, J. Geophys. 
Res., 89, 2169-2179, 1984. 

Voigt, G.-H., Magnetospheric equilibrium configurations and slow adiabatic convection, in Solar 
Wind-Magnetosphere Coupling, edited by Y. Kamide and J. A. Slavin, pp. 233-273, Terra 
Scientific Publishing Co., Tokyo, 1986. 



